module formula
	implicit none

contains
	subroutine Newton_Raphson(eps, x, answer, error,count,max)
		implicit none
		real*8 eps, x, x0, x1, x2,error,answer
		integer*2 count,max
		error = 1;
		count=0
        write(10,*) count,x,error
		do while (error > eps .AND. count < max)
			count = count+1
			x0 = x	
			call f2a(x,x1)
			call f_a2a(x,x2)
			x = x0 - (x1/x2)
			error = abs(x-x0) 
			write(10,*) count,x,error
		end do
		answer = x
	end subroutine Newton_Raphson
    
	subroutine Lelaran_NR(eps, x, y,count)
		implicit none
		real*8 eps, x, y, x0, y0,error_x,error_y
		integer*2 count
		error_x = 1
		error_y = 1
		count=0
		write (10,*) count,x,y,error_x,error_y
		do while (error_x > eps .OR. error_y > eps)
			count = count+1
			x0 = x
			y0 = y	
			call determinanJacobi(x0,y0,x)
			call determinanJacobi(x0,y0,y)
			x = x0 - x
			y = y0 + y
			error_x = abs(x-x0)
			error_y = abs(y-y0) 
			write(10,*) count,x,y,error_x,error_y
		end do
	end subroutine Lelaran_NR

	subroutine jacobi(eps, x, y,count,t)
	implicit none
	real*8 eps, x, y, x0, y0,error_x,error_y
	integer*2 count
    integer t
	error_x = 1;
	error_y = 1;
	count=0
	write(10,*) count,x,y,error_x,error_y
	do while ((error_x > eps .OR. error_y > eps) .AND. count < 100)
		x0 = x
		y0 = y
		count = count+1
        if(t .EQ. 1) then
			call fx(x0,y0,x)
			call fy1(x,y0,y)
		else if (t .EQ. 2) then
			call fx(x0,y0,x)
			call fy2(x,y0,y)
        else if (t .EQ. 3) then
			call fx(x0,y0,x)
			call fy3(x,y0,y)
        end if
		error_x = abs(x-x0)
		error_y = abs(y-y0) 
		write(10,*) count,x,y,error_x,error_y
	end do
	end subroutine jacobi

    subroutine secant (eps,a,b,hasil,count,i)
    implicit none    
    real*8 eps,a,b,hasil,x0,x1,y0,y1,y2,c
    integer*2 count
	integer i
    x0 = a
    x1 = b
    count = 0
    hasil = 0
    call fungsi (x0,i,y0)
    call fungsi (x1,i,y1)
    c = ((y0*x1)-(y1*x0))/(y0-y1)
    do while ((abs(c-x0) > eps) .AND. (count < 100))
        count = count + 1
        call fungsi (x0,i,y0)
        call fungsi (x1,i,y1)               
        if ((y0-y1) .EQ. 0.0D0) then
        write (10,*) 'Fail'
        exit
        end if       
        c = ((y0*x1)-(y1*x0))/(y0-y1)
        call fungsi(c,i,y2)       
        if (abs(y2) < eps) then
        hasil = c
        exit
        end if
		write (10,*)count,x0,x1
        x0 = x1
        x1 = c
    end do
    end subroutine secant

    subroutine bagidua (eps,a,b,jawaban,count,i)
    implicit none 
    real*8 a,b,c,jawaban,eps,deltaN,temp1,temp2,y1,y2
    integer*2 count
    integer i
    temp1 = a
    temp2 = b
    count = 0
    y1 = 0
    y2 = 0
    deltaN = abs (temp1-temp2)
    c = (temp1+temp2)/2
       
    do while ((deltaN > eps) .AND. (count < 100))
    count = count + 1
    c = (temp1+temp2)/2
    write (10,*) count,"     ",temp1,"     ",temp2,"     ",deltaN 
    
    jawaban = c
    call fungsi(temp1,i,y1)
    call fungsi(c,i,y2)
    if ((y1*y2) < 0) then
        temp2 = c
        deltaN = abs (temp1-temp2)
    else 
        temp1 = c
        deltaN = abs (temp1-temp2)
    end if    
    end do
    jawaban = c
    end subroutine bagidua
    
    subroutine regulafalsi (eps,a,b,jawaban,iter,i)
    implicit none
    real*8 a,b,eps,jawaban,c,temp1,temp2,y1,y2,y3
    integer*2 iter
    integer i
    iter = 0
    temp1 = a
    temp2 = b
    y1 = 0
    y2 = 0
    call fungsi (temp2,i,y2)
    call fungsi (temp1,i,y1)
    if (y1*y2>0) then
    write (*,*) 'Nilai tidak valid'
    iter = 101
    end if
    do while ((abs (temp1-temp2) > eps) .AND. (iter < 100))
        iter = iter +1
        call fungsi(temp2,i,y2)
        call fungsi(temp1,i,y1)
        c = temp2-( (y2*(temp2-temp1))/(y2-y1) )
        call fungsi(c,i,y3)
        write (10,*) iter,"     ",temp1,"     ",temp2,"     ",y3
        if(abs(y3) < 1.E-6) then
            temp1=c
            temp2=c
            jawaban = c
            exit  
        else if ((y1*y3) < 0) then
                temp2=c
        else
                temp1=c
        end if        
    end do
    end subroutine regulafalsi
    
	subroutine ex(x,y)
		implicit none
		real*8 x,y
		y = 1 + x + x**2/2 + x**3/6 + x**4/24
	end subroutine ex

	subroutine fx(x,y,z)
	implicit none
	real*8 x,y,z
	z = sqrt((2*y**2+1)/3)
    !z = sqrt(10-x*y)
	end subroutine fx

	subroutine fy1(x,y,z)
	implicit none	
	real*8 x,y,z
	z = (8+2*x-x**2-y**2)/2
    !z = sqrt((57-y)/(3*x))
	end subroutine fy1

	subroutine fy2(x,y,z)
	implicit none	
	real*8 x,y,z
	z = sqrt(8+2*x-x**2-2*y)
	end subroutine fy2
    
    subroutine fy3(x,y,z)
	implicit none	
	real*8 x,y,z
	z = (8+2*x-x**2)/(y+2)
	end subroutine fy3
    
subroutine f(x,z)
implicit none
real*8 x,y,z
call ex(x,y)
z = y - 5*x**2
end subroutine f

subroutine f_a(x,z)
implicit none
real*8 x,y,z
call ex(x,y)
z = y - 10*x
end subroutine f_a

	subroutine f2a(x,z)
		implicit none
		real*8 a,b,c,d,x,z
        !a = -2*x
        !b = -0.1*x
		!call ex(a,c)
        !call ex(b,d)        
		z = 80*exp(-2*x)+20*exp(-0.1*x)-10
	end subroutine f2a

	subroutine f_a2a(x,z)
		implicit none
		real*8 a,b,c,d,x,z
        !a = -2*x
        !b = -0.1*x
		!call ex(a,c)
        !call ex(b,d)        
		z = -160*exp(-2*x)-2*exp(-0.1*x)
	end subroutine f_a2a

subroutine determinanJacobi(x,y,answer)
implicit none
real*8 x,y,answer,dux,dvx,duy,dvy,u,v
call fdux(x,y,dux)
call fdvx(x,y,dvx)
call fduy(x,y,duy)
call fdvy(x,y,dvy)
call fu(x,y,u)
call fv(x,y,v)
if (answer .EQ. x) then
	answer =(((u*dvy)+(v*duy))/((dux*dvy)-(duy*dvx)))
else
	answer =(((u*dvx)-(v*dux))/((dux*dvy)-(duy*dvx)))
end if
end subroutine determinanJacobi

subroutine fu(x,y,z)
implicit none
real*8 x,y,z
z = 3*x**2 - 2*y**2 - 1
end subroutine fu

subroutine fv(x,y,z)
implicit none
real*8 x,y,z
z = x**2 -2*x + y**2 +2*y -8
end subroutine fv

subroutine fdux(x,y,z)
implicit none
real*8 x,y,z
z = 6*x
end subroutine fdux

subroutine fdvx(x,y,z)
implicit none
real*8 x,y,z
z = 2*x-2
end subroutine fdvx

subroutine fduy(x,y,z)
implicit none
real*8 x,y,z
z = -4*y
end subroutine fduy

subroutine fdvy(x,y,z)
implicit none
real*8 x,y,z
z = 2*y+2
end subroutine fdvy

subroutine fungsi(x,i,y)
INCLUDE <EXCEPTN.INS>
real*8 x,y
integer i,OldHandler, AllHandler
OldHandler = TRAP_EXCEPTION@(ALL_FLOATING_POINT,AllHandler)

if (i .EQ. 1) then
y = (x*((2.1-(0.5*x))**0.5))/((1-x)*(1.1-(0.5*x))**0.5) - 3.69
else if (i. EQ. 2) then
y = tan(x)-x+1
else if (i .EQ. 3) then
y = 0.5*exp(x/3) - sin(x)
else if (i .EQ. 4) then
y = ((x+1)/(x*(1+(0.01*x)))) - log((1+(x*0.01)) / (0.01*x))
else if (i .EQ. 5) then
y = ((x+1)/(x*(1+(0.005*x)))) - log((1+(x*0.005)) / (0.005*x))
else if (i .EQ. 6) then
y = ((x+1)/(x*(1+(0.001*x)))) - log((1+(x*0.001)) / (0.001*x))
else if (i .EQ. 7) then
y = (1.4E-5*(x**1.5)) + (1.15E-5*(x**2))
end if
end subroutine fungsi


!real*8 function f1a(x)
!implicit none
!real*8 x
!f1a = (x*((2.1-(0.5*x))**0.5))/((1-x)*(1.1-(0.5*x))**0.5) - 3.69
!end function f1a

!real*8 function f1b(x)
!implicit none 
!real*8 x
!f1b = tan(x)-x+1
!end function f1b

!real*8 function f1c(x)
!implicit none
!real*8 x
!f1c = 0.5*exp(x/3) - sin(x)
!end function f1c

!real*8 function f3a(x)
!implicit none
!real*8 x
!f3a = ((x+1)/(x*(1+(0.01*x)))) - log((1+(x*0.01)) / (0.01*x))
!end function f3a

!real*8 function f3b(x)
!implicit none
!real*8 x
!f3b = ((x+1)/(x*(1+(0.005*x)))) - log((1+(x*0.005)) / (0.005*x))
!end function f3b

!real*8 function f3c(x)
!implicit none
!real*8 x
!f3c = ((x+1)/(x*(1+(0.001*x)))) - log((1+(x*0.001)) / (0.001*x))
!end function f3c

!real*8 function f4(x)
!implicit none
!real*8 x
!f4 = (1.4E-5*(x**15)) + (1.15E-5*(x**2))
!end function f4

integer function AllHandler()
    INCLUDE <EXCEPTN.INS>
    write (10,*)'Solution not found'
    AllHandler = NONCONTINUABLE_EXCEPTION
end function AllHandler

end formula 